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On the spatial distribution of thermal energy in equilibrium 
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The equipartition theorem states that in equilibrium thermal energy is equally distributed among 
uncoupled degrees of freedom which appear quadratically in the system’s Hamiltonian. However, 
for spatially coupled degrees of freedom — such as interacting particles — one may speculate that 
the spatial distribution of thermal energy may differ from the value predicted by equipartition, 
possibly quite substantially in strongly inhomogeneous/disordered systems. Here we show that for 
systems undergoing simple Gaussian fluctuations around an equilibrium state, the spatial distribu¬ 
tion is universally bounded from above by ^ksT. We further show that in one-dimensional systems 
with short-range interactions, the thermal energy is equally partitioned even for coupled degrees of 
freedom in the thermodynamic limit and that in higher dimensions non-trivial spatial distributions 
emerge. Some implications are discussed. 


Equilibrium thermal fluctuations play a key role in 
physics, chemistry and biology, and the framework that 
captures their properties — statistical thermodynamics 
— is a central branch of physics. One of the renowned 
results obtained in this field is the equipartition theorem 
[1], which in its simplest form states that the total ther¬ 
mal energy of the system is equally distributed among its 
uncoupled degrees of freedom (DOF). In addition, each 
uncoupled DOF appearing quadratically in the Hamil¬ 
tonian has on average an energy of ^ksT, where fcs is 
the Boltzmann constant and T is the absolute tempera¬ 
ture [1]. 

The equipartition theorem holds only for uncoupled 
DOF, and strictly speaking does not state anything about 
the energy of coupled DOF. Quadratic Hamiltonians can 
always be decomposed into a set of uncoupled DOF (“mu¬ 
tually orthogonal normal modes”), and the theorem ap¬ 
plies to them. It is intriguing, though, to ask in all gen¬ 
erality whether the average potential energy of coupled 
DOF can significantly deviate from the value predicted 
by equipartition. Put simply, we ask: what can be said in 
general about the spatial distribution of thermal energy? 
To the best of our knowledge, despite its basic nature 
and generality, this question has not been posed in the 
literature, nor answered. 

One may speculate that a localized enhancement or in¬ 
hibition of thermal energy may have some effect on vari¬ 
ous local processes. For instance, if local thermal fluctu¬ 
ations activate chemical reactions, these might be facili¬ 
tated or hindered in the presence of enhanced or reduced 
fluctuations. Other processes which might be affected are 
structural changes, such as severing of biopolymers [2, 3]. 
Related effects may also be observed in elastic-network 
models of protein folding, where local (nearest-neighbor) 
fluctuations are assumed to dictate bond rupture [4]. We 
note that while much attention has been devoted to the 
effect of disorder on non-equilibrium transport proper¬ 
ties, e.g. [5-7], to the best of our knowledge none of the 
previous works explored the equilibrium spatial distribu¬ 
tion of thermal energy, as we intend to do here. 

As a prelude, we begin by solving a very simple prob¬ 
lem, depicted in Fig. I. Consider a system of two par¬ 
ticles interacting via linear springs with each other and 



FIG. 1: A model system of two masses and three linear 
springs, connected in series between fixed walls. Xi measures 
the deviation of the i-th particle from its equilibrium position. 

with bounding walls. Its potential energy is 

U = = \kix\ + \k2{x2 - Xif P \kzxl, ( 1 ) 

where Cq is the energy of the a-th spring and the xfs 
denote the deviation of the z-th particle from its equilib¬ 
rium position (when possible, we adopt the convention 
that Latin indices denote DOF, while Greek indices de¬ 
note interactions). In Eq. (I), we assume that the rest 
length of the entire chain is identical to the distance be¬ 
tween the bounding walls. As is typical in such systems, 
the kinetic energy, 'Yhi ^triixf, is a sum of quadratic un¬ 
coupled terms. Thus, diagonalizing the kinetic contribu¬ 
tion is trivial, and hence in what follows we disregard the 
kinetic energy of the system. 

What is the average energy stored in the a-th spring? 
In particular, can it significantly deviate from the value 
predicted by equipartition? For such a simple system 
the answer is readily calculable through the correlations 
between the xfs. For example, (£ 2 ) is given by 

(€ 2 ) = (^2 - = y ( (xl) + (xl) - 2 (X1X2) ) , 

where (•) denotes thermal averaging. Since the energy is 
quadratic, the correlation matrix C is given in terms of 
the Hamiltonian H by [I] 

Cij = {x,x,) = kBT{H-^),, , (2) 

where the the Hamiltonian (or the Hessian) is defined as 

j j 

Hij = g" . With these formulae, an explicit calculation 
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yields [23] 


i) = ^ksT 
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( 3 ) 


A few insights can be gained from this very simple 
example. First, a clear cut answer is given to the ques¬ 
tion presented above: for a general choice of the k^s, 
the thermally averaged energy of a given spring may dif¬ 
fer from the value predicted by equal partition. Second, 
since the system consists of two DOF (i.e. the Hamilto¬ 
nian has two normal modes), each contributes ^ksT to 
the total energy, and thus J2ai^a) — ksT, as expected. 
That is, the spatial average of the energy agrees, by 
construction, with the equipartition theorem, and reads 
{N -\- = ^ksT (with N = 2). Third, we 

note that (cq,) is bounded between 0 and which 

means that inhomogeneity in the fc^’s might either in¬ 
crease or decrease it relative to the spatially average value 
of ^ksT, depending on the inhomogeneity. Lastly, we 
note that for a homogeneous system, i.e. ka = k, the en¬ 
ergy is obviously equally partitioned among the springs 
and equals ^ksT. 

These results might appear somewhat restricted as 
they involve a system with a small number of DOF, 
N = 2, and involve specific boundary conditions which 
may play a non-trivial role, especially in a small sys¬ 
tem. Consequently, we next aim at understanding how 
the spatial distribution of thermal energy depends on the 
number of DOF, N, and on the boundary conditions. 

We first consider a system of N DOF x = {xi,... ,xn)'^ 
and -|- 1 springs. The potential energy is 


Af-l-l 




2ka (^Q l) 


( 4 ) 


where 1 <q;<A^-|-1, Ca is the potential energy of the a-th 
spring and the ka’s are non-negative constants. Formally, 
Eq. (4) makes reference to Xq or x^+i. These are not real 
DOF, but rather “boundary conditions” imposed by the 
fixed walls and should be taken as xq = xn+i = 0. The 
general question that we pose is: what can be stated 
about the distribution of (ec), for general fc^’s? 

The conventional procedure for addressing such a ques¬ 
tion is to diagonalize the Hamiltonian H and work in 
the basis of its normal modes. Then, (ea) can be recon¬ 
structed, at least conceptually, by summing over the con¬ 
tributions of the individual modes (e.g. Eq. (9) in [9]). 
While this generic recipe is very useful in most cases, in 
this case working with the normal modes obfuscates the 
structure of the problem, since there is no simple way to 
describe them for a general distribution of the k^s. 

What is then a useful basis to work with? To an¬ 
swer this question we first note that for any linear 
change of variables x = Ax, where A is an invertible 
N X N matrix, the modified Hamiltonian takes the form 
jH = A~'^HA~^, where A~'^ stands for {A~^)'^ (note 
that this is not a similarity transformation, as A is not 
orthogonal). Straightforward matrix manipulations show 


that even for the non-orthogonal variables x the correla¬ 
tion matrix is given by the inverse of the relevant Hamil¬ 
tonian, i.e. C = (xx^'^ = kBTH~^. The main insight 
gained from this brief discussion is that one should not be 
constrained to using an orthogonal basis in transforming 
the Hamiltonian into a desired form. This turns out to 
be important for solving the problem at hand. 

Following this insight, we look for new variables x = 
A X such that H will be, loosely speaking, “almost” di¬ 
agonal. A clue for finding a useful basis is obtained by in¬ 
specting Eq. (4), which is already written in an “almost” 
diagonal form, if we identify the new variables simply 
as Xi = Xi — Xi-i- That is, we take the combinations 
that make up the interactions as the new variables. This 
defines the transformation matrix Aij = dij — Sij^i. 

Under this choice of non-orthogonal variables, almost 
all of the spring energies in Eq. (4) become Cq, = ^kaX^- 
The last relation, however, is valid only for 1 < a < A^. 
Clearly, a one-to-one correspondence between the DOF 
and the interactions/springs is impossible since the num¬ 
ber of springs exceeds N. Indeed, an explicit calculation 
shows that in terms of the new variables the energy is 
not strictly decoupled, but only almost 


U = Ea=l 5*0^0 + lkN+1 (il H-h XNf 

Hij = ktSij + kN+i , or, H = K + kN+ibb^ , 


( 5 ) 


where K = dia,g{ki,..., kN) and b is a vector of I’s. 

Equation (5) is very useful since the inverse H~^ 
is readily calculated using the Sherman-Morrison for¬ 
mula [10], which can be expressed in the form 


(H + kvv^) ^ =H-^ 


H-^vv'^ i/ 1 
k~^ + v'^H~^v 


( 6 ) 


and is valid whenever both H and H + kvv"^ are in¬ 
vertible (here v is a vector and fc is a scalar). Applying 
this formula to Eq. (5), we obtain the generalization of 
Eq. (3) to any N 


(Ca) = IksT 


1-iN + l) 
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( 7 ) 


where ((fc ^)) = (A^-|-1) is the quenched aver¬ 

age. 

Two features of this result will prove important. First, 
it is seen that (ca) < ^ksT regardless of the choice of ka’s 
(moreover, the order in which the ka ’s are distributed in 
space makes no difference, (eq) depends only on k~^ and 
the average ((fc“^))). Second, it is evident that (£„) tends 
towards ^ksT in the thermodynamic limit N —>■ oo, as 
long as ka^/{{k~^)) does not increase with N. 

Our next task is to generalize Eq. (7). It will be shown 
that these two features are general for a wide class of 
physical systems, namely systems with local interactions 
undergoing Gaussian fluctuations around a stress-free 
equilibrium. As we will show, for these systems ^ksT is 
a strict upper bound for (cq). This result is entirely gen¬ 
eral, independent of dimensionality or interaction range. 
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Second, the fact that (cq) = ^ksT plus a negative cor¬ 
rection of order N~^, which depends on the inhomogene¬ 
ity, is the general rule for one-dimensional systems with 
short-range interactions. In particular, in such systems 
the distribution of thermal energy becomes spatially con¬ 
stant in the thermodynamic limit. These are two main 
results of this work. 

Consider a system with N DOF and total energy U = 
where n is the number of interactions. The 
most general expansion of £„ in the DOF reads 

= + + . ( 8 ) 

The linear term vanishes under thermal averaging as long 
as anharmonic contributions to the energy are neglected, 
and hence is omitted hereafter (note that = 

dU/dxi vanishes due to global equilibrium). The only 
assumption we adopt is that Xj can be writ¬ 
ten as Xi)^, where S is an nx matrix which 

describes the interactions in the system (for example, B 
can be easily read off Eq. (4)). In this case, the Hamil¬ 
tonian is given by H — B^B. 

This is a generic form of local interaction energies for a 
broad class of physical systems: (i) Discrete field theories, 
or discrete approximations to continuum field theories, 
where the energy density takes the form [£(/)]^, with 
some spatial linear differential operator C and field /. 
Relevant examples include — among many others — the 
Euler-Bernoulli theory of elastic beams [11], the Foppl 
von-Karman theory of thin sheets [11] and the Helfrich 
theory of membrane elasticity [12]. (ii) Systems of dis¬ 
crete particles interacting via a radially-symmetric pair¬ 
wise potential, where at equilibrium all particle pairs are 
at a stress-free configuration. Relevant examples include 
glassy systems near jamming [13] and elastic networks [4]. 

Using Eq. (2), (eq) can be readily expressed in terms 
of the interaction matrix B, as 

(ca) = P = B{B^By^B^. (9) 

P is an orthogonal projection operator [14], since P^ = P 
(this holds even when B^B is not invertible [15]). A gen¬ 
eral property of such operators is that all their elements 
are smaller than unity in absolute value. We thus prove 
a central result of this work, i.e. that 

(Ca) = IksT Paa < ^bT . ( 10 ) 

The rank of P, which equals that of H and B, carries 
important information. For example, since P is a projec¬ 
tion operator that works in a space of dimension n, in case 
that rank P = n we can immediately conclude that P is 
the identity and thus (e^,) = ^bT identically. This hap¬ 
pens whenever the rows of the interaction matrix B are 
linearly independent (and in particular n<N). The case 
n = N, for nearest-neighbor interactions, corresponds to 
isostatic systems, i.e. systems where the number of con¬ 
straints (interactions) equals the number of DOF [13]. 


Equation (10) puts strict bounds on the possible values 
of (cq,), but much more can be said about the behavior 
within these bounds. Specifically, we can derive the ana¬ 
log of Eq. (7) in the general case of one-dimensional sys¬ 
tems with short range interactions. The detailed deriva¬ 
tion can be found in [15], and it follows verbatim the 
structure of the derivation of Eq. (7), as outlined here: 
For one-dimensional systems, where each DOF interacts 
with its m nearest neighbors the number of interactions 
generally exceeds the number of DOF by m. A non- 
orthogonal transformation is used to bring the Hamil¬ 
tonian to the form H = K -\- kabab^, where K is 
diagonal and the second term is a sum over the m “ex¬ 
cess” interactions. Then, the Sherman-Morrison formula 
is iteratively applied m times to calculate the inverse. Be¬ 
cause of the short-range nature of the interactions, the 
magnitude of the non-diagonal correction to H~^ is of 
order m/N, and vanishes in the thermodynamic limit. 
Thus, (ea) = yBT + 0{N-^). 

This is another central result of this work: in one¬ 
dimensional systems with short range interactions, the 
spatial distribution of thermal energy becomes essentially 
flat in the thermodynamic limit. The crux of the argu¬ 
ment lies in the fact that the number of interactions does 
not greatly exceeds the number of DOF, and that the 
ratio between them approaches unity in the the thermo¬ 
dynamic limit. 

Returning to the problem of the spring-mass chain, and 
having solved the problem of the A^-dependence for fixed 
boundary conditions, we now turn to explore the effect 
of boundary conditions for a fixed N. For instance, semi¬ 
fixed boundary conditions can be obtained by removing 
the constraint of one of the walls, by setting, say, ki=0. 
In doing so we obtain n = N independent interactions and 
thus (cq) = is identically constant for any N. Fully 

free boundary conditions are obtained by setting both 
ki and k^+i to zero, and give rise to a single Goldstone 
mode (uniform translation). In this case we have n = N — 
I independent interactions and again (cq.) is identically 
constant (clearly, both these results can also be obtained 
by properly taking limits of Eq. (7)). This also shows 
that in one-dimensional systems the effect of boundary 
conditions in non-local, i.e. every spring in the system is 
affected by the bounding walls. 

The general approach discussed above can be applied 
to different types of interactions. For example, bend¬ 
ing fluctuations are described by local interactions of the 
form Cq = {xa-i — ‘2.Xa+Xa+i)^^ where the Kq,’s are 
the local bending rigidities. Identical arguments show 
that in a chain with free boundary conditions, the spa¬ 
tial distribution of bending fluctuational energy is exactly 
constant, regardless of the choice of the Kq’s. 

This result offers the first application of the theoreti¬ 
cal development described in this paper, as it seems to 
refute a recently conjectured mechanism for severing of 
actin filaments, one of the most important and ubiq¬ 
uitous biopolymers in eukaryotic cells [16]. In [17-19], 
it has been hypothesized that thermal energy may be 
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concentrated at the boundaries between relatively softer 
and stiffer regions of the biopolymer (softening is in¬ 
duced by a different molecule, cofilin, which partially 
binds actin [17, 20]), and that the excess thermal energy 
is responsible for the experimentally observed preferen¬ 
tial severing near these boundaries. Our result shows, 
at least within the framework of a discrete description 
of quadratic bending fluctuations, that no such energy 
concentration takes place. 

What happens in higher dimensions? As the crux of 
the argument lies in the relative number of DOF and in¬ 
teractions, dimensionality appears to be crucial. In fact, 
in dimensions higher than one the argument seems to fail 
qualitatively as generally there are significantly more in¬ 
teractions than DOF. In this case, we also expect local 
topological variations, bond strength disorder, defects, 
holes, free boundaries and the like to play a role. 

To see this, consider a hexagonal portion of a two- 
dimensional triangular lattice of particles (which 
amounts to N DOF), interacting via linear springs, as 
shown in the inset Fig. 2. Clearly, in the limit of large 
systems the number of springs n approaches |fV. When 
all of the springs are identical, i.e. with no bond strength 
disorder, we expect the average energy of a spring far 
from the free boundary to approach {ea) ^ ^ksTx ^ ^ 
(this was verified by an explicit calculation). There 
is no reason, however, to expect the thermal energy to 
be spatially uniform in the presence of inhomogeneities, 
either in the lattice topology or in the bond strength. 

To test this, we considered the lattice in the inset of 
Fig. 2 with bond strength disorder, where the k^’s are 
normally distributed [24]. (ca) is plotted vs. ka in the 
main panel, where the average energy of bulk springs in a 
homogeneous system, ^ksT, is shown as well. Bulk and 
boundary springs are distinguished. Several key observa¬ 
tions can be made: (i) Unlike in one-dimensional systems, 
thermal energy spans the whole interval between 0 and 
both below and above the homogeneous system 
bulk value, ^ksT. (ii) (ca) appears to vary systemati¬ 
cally with the local spring strength ka- (in) Boundary 
springs have higher energy than bulk springs. This is a 
purely topological effect in which boundary springs have 
less neighbors than bulk springs, an effect that persists 
near free boundaries in fully ordered systems. In gen¬ 
eral, disordered systems (e.g. glassy ones [13]) feature 
also bulk topological disorder. 

In summary, in this work we posed a basic question 
in statistical physics: What is the spatial distribution 
of thermal energy in equilibrium? We showed that un¬ 
der the stated conditions it is strictly bounded between 
0 and and that for one-dimensional systems with 

short-range interactions, the spatial distribution of ther¬ 
mal energy becomes essentially flat in thermodynamic 
limit, even for highly disordered systems. The crux of 
the derivation lies in the fact that in one-dimensional 
systems the number of interactions is the same as the 
number of DOF, up to an additive constant which is neg¬ 
ligible in the thermodynamic limit. In higher dimensions 



FIG. 2: (Color online) The average thermal energy {sa.} (in 
units of ksT) as a function of ka/ {{k}), for a hexagonal por¬ 
tion of a two-dimensional triangular lattice. The data are 
partitioned into bulk (blue) and boundary springs (yellow), 
cf. inset. Spring constants are distributed normally with 
mean 1 and variance 0.3. Each side of the hexagon consists 
of 20 springs, which means N = 2522 (a smaller system, with 
3 springs and N = 74, is shown in the inset for illustration). 
No quantitative change was observed with increasing N. The 
points show data from 30 realizations and the solid lines are 
guides to the eye. The dashed line shows the asymptotic 
value |, corresponding to (ea) in the bulk of a homogeneous 
triangular lattice in the thermodynamic limit. 


this does not hold, as was explicitly demonstrated in a 
specific example. Systematically unraveling the relations 
between the spatial energy distribution and dimensional¬ 
ity, the system’s geometry and the form of disorder is a 
theoretical challenge for future work. 

The most outstanding question that emerges from this 
work is what the influence of the spatial distribution of 
thermal energy on various physical processes and quanti¬ 
ties might be. If local energy fluctuations can affect local 
process, as was suggested — for example — in the con¬ 
text of bond rupture in elastic-network models of protein 
folding [4], then one can imagine the possibility of tai¬ 
loring high dimensional systems in order to enhance or 
reduce thermal fluctuations in defined locations, to con¬ 
trol local processes of interest. Another important future 
direction would be to explore the roles played by stresses 
(both internal and external). 
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This short document is meant to provide additional 
technical details in relation to results appearing in the 
manuscript. 


I. GENERAL FORMALISM 


The general description of the problem addressed in 
the manuscript is as follows. Consider a system of N de¬ 
grees of freedom x = {xi,... and n local quadratic 

interactions between the DOF, 

e„ = Y Da^x^ =^{da-xf , l<i<n 

(SI) 

where Z3 is a dimensionless nxN matrix whose a-th row 
is da, ka are non-negative numbers with dimensions of 
energy and x is dimensionless. The non-zero elements 
of D are assumed to be of order unity. Note that this 
notation and units are slightly different compared to the 
manuscript as our aim here is to remain as general as 
possible. The results reported in the manuscript are re¬ 
covered when we identify ba = \/k^da and B — Vk^D, 
where = ka 5a0 is an n x n diagonal matrix. Fol¬ 
lowing the discussion below Eq. (8) in the main text, the 
linear term of Eq. (SI) is omitted. 

The total energy is obtained by summing over Cq, and 
the Hamiltonian can be expressed in terms of D as 

n 

H = D^K^^^D = '^ka dadl . (S2) 

CK —1 

As described in the main text, (cq) is determined by the 
diagonal elements of the projection operator 

P = , (S3) 

where f stands for the Moore-Penrose pseudo-inverse, see 
Sect. II for definitions and details. For an invertible ma¬ 
trix, f is simply the inverse. 

In order to proceed, there are two distinct cases to 
examine: 


(i) {da} are linearly independent (and hence n<N). 

(ii) {da} are linearly dependent. 

Case (i) is relatively simple. In this case the rank of 
B is n, and so is the rank of P. Since P is a projection 
operator in an n-dimensional space, we conclude that P 
must be the identity operator, and thus (cq) = ^ksT. 

It might be useful to explicitly construct the trans¬ 
formation matrix A that shows how this happens. By 
the rank-nullity theorem we know that there are N — n 


vectors g„+i, ■ ■ ■ ,gN (“Goldstone modes”) which satisfy 
Dga=0 and do not have an energetic cost. The n rows 
of D together with the N — n Goldstone modes form a 
basis, i.e. the new variables take the form x = Ax with 


/- di -\ 


- gn+l - 

V- gN -/ 


(S4) 


In these terms, the energy takes the form 



(S5) 

With this at hand it is immediate to see that 

(eo) = Y (ia) = y^bT = \kBT, (S6) 

since 


H^ = 




(S7) 


where = ka ^ ^a/s- 

We thus conclude case (i) with the following out¬ 
come: whenever the interactions (either ba or da) are 
linearly independent, all of the local energies are equal, 
(ca) = ^ksT, irrespective of any other property of the 
interactions. The total energy of the system is then 
^nksT, which amounts to assigning ^ksT to every non- 
Goldstone mode. 

Case (ii) is more involved in the sense that it is difficult 
to make further progress without substantially constrain¬ 
ing the structure of the interactions. As will be shown, 
in order for the result to remain valid (at least in the 
thermodynamic limit) the number of interactions must 
not greatly exceed the number of DOF. This is always 
the case in one-dimensional systems with short range in¬ 
teractions, as we now show. 

In the following we consider a one-dimensional systems 
with local interactions of the form 



1 < a < N + m (S8) 


where to G N is the interaction length. Since we are in 
one-dimension, in this case we have n = N + m. We 
remind the reader that this is a shorthand writing where 
the summand might formally make reference to Xi with 
i < 0 or i> N. These are to be understood as externally 
specified “boundary conditions”, and should be taken as 
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0. For simplicity, we assume here that H is invertible, 
but this assumption is not essential (see Sect. II). 

For these interactions, it is guaranteed that di,..., djq 
are linearly independent, and we can thus dehne 

Xa = da ■ X i < cx < N , (S9) 

or equivalently x = DqX, where Dq is the invertible 
matrix composed of the first N rows oi D. In terms of 
the new variables, the Hamiltonian takes the form 


n 

H = K^^'>+ kaD^^dadlD^^ 

ol—N-\-1 

n 

= FS:W+ ^ kadadl , 


(SIO) 


with the notation da = Dg da- Note that here = 
ka^a 0 is a NxN matrix. 

It is seen that ff is a diagonal matrix, supplemented 
with m non-diagonal corrections (recall that n = N + m). 
Since in the absence of these corrections the result (cq) = 
^ksT follows immediately (for 1 < a < N), it is left 
only to verify that they do not significantly affect H~^. 
To this end, we iteratively make use of the Sherman- 
Morrison formula [S5] which can be put in the form 

{H + kvv^)~^ =H-^- ^ ' , (SIl) 

'' ' k-^+v^H-^v ^ ’ 

and is valid whenever both H and H + k vv^ are invert¬ 
ible (here u is a vector and A: is a scalar). To calculate the 
inverse of H, as given in Eq. (SIO), we iteratively apply 
the formula m times, as follows. We define Hq = 
and 


Ha — Ha —1 T k]\[-\.a 1 ^ ^ ^ ■ (SI2) 


Each iteration, is calculated in terms of Ha and 

da+i using Eq. (SlI). 

We now turn to estimate the corrections to the inverse, 
i.e. the last term in Eq. (SIl) in each iteration. Since Dq 
is a sparse matrix, having non-zero elements only on m 
diagonals, is generically a dense matrix. In fact, in 
the relatively simple case discussed here, a closed form 
recursive expression can be obtained, reading (£)g"^)^^. = 
/(* - j), with 


z < 0 

fiz) = <1 Co z = 0 

-Co ^ 1 fiz - k)ck z > 0 


(S13) 


as can be readily verihed by explicit calculation. Thus, 
the number of non-zero elements of da scales with N (as 
opposed to da, for which this number is exactly m). Ex¬ 
plicitly, the (a, /?) element of the non-diagonal correction 
in the first iteration takes the form 

ka k 1^ da d^ 


AT-1-1 


■ E,"Li k. 


O' d^d'y 


where da is shorthand for the a-th component of d]\[^i, 
which is of order unity. Note that no summation over 
a, (3 is implied. The above expression shows that the de¬ 
nominator of this correction contains 7V-I-I terms of order 
{{k~^)) and hence it scales as N ((^k~^)). The numerator, 
obviously, has no N dependence. Consequently, we have 

+ ^+0{N-‘^) (S14) 

where A is a matrix whose elements are of order 
ka^kj^/ (^k~^'f). This argument can be similarly applied 
in all m iterations, as long as m N, which is always 
valid in the thermodynamic limit. We can thus conclude 
that the correlation matrix is 

C = kBTH-^ = kBTK-^ z,o{N-^) , (S15) 

and that the average value of the local energies is 

(e.) = ^kBT (h-^) = \kBT + O(iV-i) . (S16) 

Z \ / act Z 

II. SUM OVER NON-GOLDSTONE MODES 

In the manuscript, and in the previous section, it 
was stated that calculating correlations by summing over 
non-Goldstone modes is equivalent to using the Moore- 
Penrose inverse [SI-S3] of H. This follows trivially from 
the definition of the latter. 

The Moore-Penrose pseudo-inverse is a generalization 
of the usual inverse of a matrix. Eor a symmetric square 
NxN matrix A, the Moore-Penrose inverse — denoted 
by At — can be explicitly calculated in terms of the 
eigenvectors of A [S2, S3] according to 

4 = 5: , (S17) 

where Q diagonalizes A, in the sense that QHQ^ is 
diagonal and the i-th row of Q is an eigenvector with an 
associated eigenvalue A^. 

This formula might appear more transparent in bra-ket 
notation, in which it reads 

At = ^ MM , (S18) 

where the |(j'a)’s form an orthonormal basis (i.e. the rows 
of Q). It is easy to verify by explicit calculation that this 
formula satisfies the four equations defining the pseudo¬ 
inverse [S2, page 40]. In this writing it is also clear that 
when A is invertible, and hence all of the A^’s are non¬ 
zero, we have At = A~^. 

Consider now a system of N degrees of freedom, x = 
{xi ,..., xn)'^, with an Hamiltonian H. We work in the 
basis of normal modes, and diagonalize H with an or¬ 
thonormal matrix Q, in the sense defined above. The 
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eigenmodes qi,..., qn are the rows of Q. The equipar- 
tition theorem [S4] says that 

{Qi ■ Qj) = > (S19) 

whenever Ai ^ 0. From this, the correlation matrix C of 
the original DOF can be calculated by summing over the 
eigenmodes 


In the case where H is invertible, we immediately iden¬ 
tify Sais/Xa with (the elements of) QH~^Q^, and thus 
Eq. (2) in the manuscript is obtained. In case some of 
the Ai’s vanish, the sum is preformed only on the non¬ 
vanishing ones and the result identifies with Eq. (SI7). 


{XiXj) — Qki 

k,l 



(S20) 


[51] R. Penrose, Math. Proc. Cambridge Philos. Soc. 51, 406 
(1955). 

[52] A. Ben-Israel and T. N. E. Greville, Generalized inverses 
(Springer-Verlag, 1974). 

[53] C. D. Meyer, Matrix analysis and applied linear algebra, 


(SIAM, 2000). 

[S4[ K. Huang, Statistical mechanics (Wiley, 1963). 

[S5] J. Sherman and W. J. Morrison, Ann. Math. Stat. 21, 
124 (1950). 



